R Notebook

This script reads in both PBMC and Liver Seurat objects, and generates select figures used in the manuscript.

Prepare objects

Load libraries

library(Seurat)

Loading required package: SeuratObject
Loading required package: sp

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

    intersect, saveRDS

Loading Seurat v5 beta version 
To maintain compatibility with previous workflows, new Seurat objects will use the previous object structure by default
To use new Seurat v5 assays: Please run: options(Seurat.object.assay.version = 'v5')

library(scClustViz)

Loading required package: shiny

library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(rcartocolor)
library(SeuratWrappers)

Read in liver map

sobj <- readRDS("~/Dropbox/Zoe/scf_version/analysis/healthy_sc/seurat_objects/dropletQC_filtered/allIntegrated_cca_kanchor5_noBiopsyHeps_dropletQCFiltered.RDS")
res <- "integrated_snn_res.1.4"
Idents(sobj) <- res
tissue <- "liver"

Read in PBMC map

load("~/Dropbox/Zoe/scf_version/analysis/healthy_sc/seurat_objects/no_dropletQC/integrated_PBMC_cca_kanchor5_scClustViz.RData")
sobj <- scSeurat
res <- "integrated_snn_res.0.6"
Idents(sobj) <- res
tissue <- "PBMC"

Visualization of metadata

UMAP with no cluster numbers

plot <- DimPlot(sobj, label = FALSE) & NoLegend()
plot
pdf(paste("./figures/", tissue, "/", tissue, "_UMAP_clusters_noLabels.pdf", sep = ""))

plot
dev.off()

png 
  2 

UMAP with cluster numbers

plot <- DimPlot(sobj, label = TRUE)
plot
pdf(paste("./figures/", tissue, "/", tissue, "_UMAP_clusters_labels.pdf", sep = ""))

plot
dev.off()

png 
  2 

Map with SCINA-generated cell-type labels

DimPlot(sobj, group.by = "scina_labels_refined", label = TRUE) & NoLegend()

Map with general cell-type labels for paper

plot <- DimPlot(sobj, group.by = "general_cell_labels",
        label = TRUE, repel = TRUE) +
  ggtitle(NULL)
plot
pdf(paste("./figures/", tissue, "/", tissue, "_UMAP_general_cell_labels.pdf", sep = ""),
    height = 7,
    width = 9)

plot
dev.off()

png 
  2 

Map grouping by general cell type labels but with no labels on plot

plot <- DimPlot(sobj, group.by = "general_cell_labels") +
  ggtitle(NULL) &
  NoLegend()
plot
pdf(paste("./figures/", tissue, "/", tissue, "_UMAP_general_cell_labels_noLabels.pdf", 
          sep = ""))

plot
dev.off()

png 
  2 

Map with original identities

plot <- DimPlot(sobj, group.by = "orig.ident",
        cols = carto_pal(length(levels(as.factor(sobj$orig.ident))), "Safe")) +
  ggtitle(NULL)
plot
pdf(paste("./figures/", tissue, "/", tissue, "_UMAP_orig_idents.pdf", 
          sep = ""))

plot
dev.off()

png 
  2 

Barplot with original identities on a cluster-level grouping:

# Meta data to plot:
df <- sobj@meta.data
# Check what column the cluster identities are in
col <- which(colnames(df) == res)
# Order clusters
df[,col] <- factor(Idents(sobj),
                             levels = c(sort(as.numeric(levels(Idents(sobj))))))
# Basic plot of clusters by replicate
ggplot(df, aes(x = get(res), fill = orig.ident)) +
  geom_bar() +
  theme(axis.text = element_text(size = 7))

# Plot as proportion or percentage of cluster
ggplot(df, aes(x = get(res), fill = orig.ident)) +
  geom_bar(position = "fill") +
  theme(axis.text = element_text(size = 7))

Barplot with original identities on a grouped by general cell labels:

df <- sobj@meta.data
plot1 <- ggplot(df, aes(x = general_cell_labels, fill = orig.ident)) +
  geom_bar() +
  scale_fill_carto_d(name = NULL, palette = "Safe") +
  theme_bw() +
  theme(axis.text = element_text(size = 8),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        axis.title.x = element_blank()) +
  ylab("Number of cells")
plot1
pdf(paste("./figures/", tissue, "/", tissue, "_barplot_orig_ident_counts.pdf", 
          sep = ""))

plot1
dev.off()

png 
  2 

# Plot as proportion or percentage of cluster
plot2 <- ggplot(df, aes(x = general_cell_labels, fill = orig.ident)) +
  geom_bar(position = "fill") +
  scale_fill_carto_d(name = NULL, palette = "Safe") +
  theme_bw() +
  theme(axis.text = element_text(size = 8),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        axis.title.x = element_blank()) +
  ylab("Proportion of cells")
plot2
pdf(paste("./figures/", tissue, "/", tissue, "_barplot_orig_ident_proportions.pdf", 
          sep = ""))

plot2
dev.off()

png 
  2 

Dotplots

Generate dotplot with specific markers

DotPlot(sobj,
        assay = "SCT",
        features = c("PTPRC", "CALCRL", "NKG7", "CD3E", "MARCO", "LYZ-1", "CD19", "MS4A1", "STAB2")
        ) +
  ggtitle("Select features for liver map")

Heatmaps

Calculate markers for general cell labels then reset resolution

Idents(sobj) <- "general_cell_labels"
sobj_markers <- RunPrestoAll(sobj, 
                             only.pos = TRUE, 
                             min.pct = 0.25, 
                             logfc.threshold = 0.25)

Calculating cluster B cells
Calculating cluster Lymphocytes 3
Calculating cluster Monocytes 1
Calculating cluster Erythrocytes
Calculating cluster Lymphocytes 1
Calculating cluster Monocytes 2
Calculating cluster Lymphocytes 2
Calculating cluster Lymphocytes 5
Calculating cluster Monocytes 3
Calculating cluster Monocytes 4
Calculating cluster Plasma B cells
Calculating cluster Unknown
Calculating cluster Mast cells
Calculating cluster Lymphocytes 4
Calculating cluster Hematopoietic stem cells
Calculating cluster Megakaryocytes

sobj_markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)

Idents(sobj) <- res

Save markers

groups = "general_cell_labels"
write.table(sobj_markers,
            file = paste("./figures/", tissue, "/",
                         tissue, "_markers_", groups, ".tsv", 
                         sep = ""),
            quote = FALSE,
            sep = "\t",
            row.names = FALSE,
            col.names = TRUE)

Generate heatmap with top 5 markers grouping by general cell types

# Remove mikado genes from marker list
sobj_markers <- sobj_markers[grep("mikado", rownames(sobj_markers), invert = TRUE),]
sobj_markers %>%
    group_by(cluster) %>%
    top_n(n = 5, wt = avg_log2FC) -> top
# If liver, select fewer cells
if (tissue == "liver") {
  cells <- sample(colnames(sobj), size = 30000)
} else if (tissue == "PBMC") {
  cells <- colnames(sobj)
}
DoHeatmap(sobj, features = top$gene, group.by = "general_cell_labels", size = 3, 
          angle = 90, cells = cells) +
  NoLegend() +
  theme(text = element_text(size = 7))

Make PDF of heatmap

groups <- "general_cell_labels"
pdf(paste("./figures/", tissue, "/", tissue, "_heatmap_", groups, ".pdf", sep = ""),
    height = 11,
    width = 7)
DoHeatmap(sobj, features = top$gene, group.by = "general_cell_labels", size = 2, 
          angle = 90, cells = cells) +
  NoLegend() +
  theme(text = element_text(size = 7))
dev.off()

null device 
          1 

Feature plots

Make specific plots with specific genes. The genes we are interested in include: PTPRC, CALCRL, NKG7, CD3E, MARCO, LYZ, CD19, MS4A1, STAB2, ALB, CD4, CD8A, CLEC4G, CD5L, C1QB, ACTA2, VWF, IGLL5, CD68. Can also plot in italics.

geneCode <- "sct_LYZ-1" # Woodchuck-specific nomenclature for this genome
gene <- "LYZ"
mapType <- "Liver"
FeaturePlot(sobj, features = geneCode) +
  ggtitle(paste(gene, "-", mapType, "map"))

FeaturePlot(sobj, features = geneCode) +
  ggtitle(bquote(~italic(.(gene))))

Another version of the feature plot that outputs genes in italics

if (tissue == "PBMC") {
 geneCodes <- c("sct_PTPRC","sct_NKG7","sct_CD14",
               "sct_CD3E","sct_MARCO","sct_LYZ-1",
               "sct_CD19","sct_MS4A1","sct_STAB2",
               "sct_CD4","sct_CD8A", "sct_XCL1;XCL2",
               "sct_CD5L","sct_C1QB", "sct_LEF1",
               "sct_ACTA2","sct_VWF","sct_IGLL5-1",
               "sct_CD68","sct_FCGR3A;FCGR3B","sct_TOP2A")
 genes <- c("PTPRC","NKG7","CD14","CD3E","MARCO","LYZ",
           "CD19","MS4A1","STAB2","CD4","CD8A", "XCL1;XCL2",
           "CD5L","C1QB","LEF1","ACTA2","VWF","IGLL5","CD68",
           "FCGR3A;FCGR3B","TOP2A") 
} else if (tissue == "liver") {
  geneCodes <- c("sct_Ptprc","sct_CALCRL","sct_NKG7",
                 "sct_CD3E","sct_MARCO","sct_LYZ-1",
                 "sct_CD19","sct_MS4A1","sct_STAB2",
                 "sct_ALB-1","sct_CD4","sct_CD8A",
                 "sct_CLEC4G","sct_CD5L","sct_C1QB",
                 "sct_ACTA2","sct_VWF","sct_IGLL5-1",
                 "sct_CD68","sct_XCL1;XCL2","sct_LEF1",
                 "sct_RSPO3","sct_MECOM")
  genes <- c("PTPRC","CALCRL","NKG7","CD3E","MARCO","LYZ",
             "CD19","MS4A1","STAB2","ALB","CD4","CD8A",
             "CLEC4G","CD5L","C1QB","ACTA2","VWF","IGLL5","CD68",
             "XCL1;XCL2","LEF1","RSPO3","MECOM") 
}
for(num in 2:length(geneCodes)) {
  plot <- FeaturePlot(sobj,
                      features = geneCodes[num]) +
    ggtitle(bquote(~italic(.(genes[num]))))
  print(plot)
  pdf(paste("./figures/", tissue, "/", tissue, "_", genes[num], "_UMAP.pdf", sep = ""))
  print(plot)
  dev.off()
}

Cell-population-specific markers

Plot marker genes for specific cell populations

Endothelial cells:

geneCodes <- c("sct_STAB2","sct_ITGA1","sct_CD55","sct_LYVE1",
               "sct_CD34","sct_VWF","sct_IFITM3;IFITM2;IFITM1","sct_RSPO3",
               "sct_MECOM","sct_Mecom","sct_CALCRL","sct_LOC114089654",
               "sct_RAMP2","sct_BST2","sct_CLEC4G","sct_CTSV;CTSL",
               "sct_STAB1","sct_PLAC9","sct_PECAM1","sct_WNT2",
               "sct_MYL6","sct_NR2F2")
genes <- c("STAB2","ITGA1","CD55","LYVE1",
           "CD34","VWF","IFITM3","RSPO3",
           "MECOM","Mecom","CALCRL","LOC114089654",
           "RAMP2","BST2","CLEC4G","CTSV;CTSL",
           "STAB1","PLAC9","PECAM1","WNT2",
           "MYL6","NR2F2") 
for(num in 1:length(geneCodes)) {
  plot <- FeaturePlot(sobj,
                      features = geneCodes[num]) +
    ggtitle(bquote(~italic(.(genes[num]))))
  print(plot)
  pdf(paste("./figures/", tissue, "/", tissue, "_", genes[num], "_UMAP.pdf", sep = ""))
  print(plot)
  dev.off()
}

Lymphocytes:

geneCodes <- c("sct_CD3D","sct_CD3E","sct_CD8A","sct_NKG7",
               "sct_KLRD1","sct_KLRB1","sct_GZMA","sct_LEF1",
               "sct_IL7R-1","sct_EOMES","sct_TIGIT","sct_XCL1;XCL2",
               "sct_TOX","sct_CTLA4","sct_CD4","sct_GZMK")
genes <- c("CD3D","CD3E","CD8A","NKG7",
           "KLRD1","KLRB1","GZMA","LEF1",
           "IL7R","EOMES","TIGIT","XCL1;XCL2",
           "TOX","CTLA4","CD4","GZMK") 
for(num in 1:length(geneCodes)) {
  plot <- FeaturePlot(sobj,
                      features = geneCodes[num]) +
    ggtitle(bquote(~italic(.(genes[num]))))
  print(plot)
  pdf(paste("./figures/", tissue, "/", tissue, "_", genes[num], "_UMAP.pdf", sep = ""))
  print(plot)
  dev.off()
}

Mesenchyme:

geneCodes <- c("sct_HHIP","sct_COL1A2","sct_COL3A1","sct_IGFBP7",
               "sct_IGFBP3","sct_DCN","sct_COL1A1","sct_SPARC",
               "sct_RBP1", "sct_CALCRL")
genes <- c("HHIP","COL1A2","COL3A1","IGFBP7",
           "IGFBP3","DCN","COL1A1","SPARC",
           "RBP1","CALCRL") 
for(num in 1:length(geneCodes)) {
  plot <- FeaturePlot(sobj,
                      features = geneCodes[num]) +
    ggtitle(bquote(~italic(.(genes[num]))))
  print(plot)
  pdf(paste("./figures/", tissue, "/", tissue, "_", genes[num], "_UMAP.pdf", sep = ""))
  print(plot)
  dev.off()
}

Look at zonation gene signatures:

cv_genes <- c("FETUB","HMGCS1","CYP2E1","GLUD1;GLUD2","CYP1A2",
              "RGN","INMT","COMT","FDPS","SPR-1")
cv_genes

 [1] "FETUB"       "HMGCS1"      "CYP2E1"      "GLUD1;GLUD2" "CYP1A2"      "RGN"        
 [7] "INMT"        "COMT"        "FDPS"        "SPR-1"      

pp_genes <- c("Saa2;Saa1-1","HAMP","APOA1","APOC2","CRYL1",
              "AMY1A;AMY1C;AMY1B;AMY2A;AMY2B","MT-ATP6","UROC1",
              "APOA2","MT-CO3")
pp_genes

 [1] "Saa2;Saa1-1"                   "HAMP"                         
 [3] "APOA1"                         "APOC2"                        
 [5] "CRYL1"                         "AMY1A;AMY1C;AMY1B;AMY2A;AMY2B"
 [7] "MT-ATP6"                       "UROC1"                        
 [9] "APOA2"                         "MT-CO3"                       

Correlation plots

Generate heatmaps comparing woodchuck clusters with various datasets

Setup for all correlations

First, set up whatever woodchuck dataset I am working with by reading in the ortholog table, choosing whether the orthologs to use are human, mouse, or woodchuck, and calculating the average expression for each cluster. The output of this section is the scaled cluster gene-expression matrix

groups <- "general_cell_labels"
# Read in ortholog table
geneNameTable <- read.table("~/Dropbox/Zoe/scf_version/make_gtf/orthofinder_sc2/homologene/collectedOrthofinderPairings.tsv",
                            sep = "\t",
                            header = TRUE)
woodchuckClusterAverages <- AverageExpression(sobj,
                                              assays = "SCT",
                                              slot = "scale.data",
                                              group.by = groups)
# Scale data
woodchuckClusterAverages$SCT <- na.omit(t(scale(t(as.matrix(woodchuckClusterAverages$SCT)))))
# Grab gene names from Seurat object
uniqueHier <- row.names(woodchuckClusterAverages$SCT)
uniqueHier <- as.data.frame(uniqueHier)
# Bind with geneNameTable to get correct order (notice uniqueHier is on left)
newNames <- dplyr::left_join(uniqueHier, geneNameTable, by = "uniqueHier")
# Get orthologs from either mouse or human
species <- "human"
if (species == "human") {
  # If human one-to-one ortho has NA, replace with mikado_final_sc2_stringent_noMito_protein column
  # This is to avoid and potential mistakes in recognizing things it shouldn't be recognizing
  newNames$speciesOneToOne <- ifelse(is.na(newNames$humanOneToOne), newNames$uniqueHier, newNames$humanOneToOne)
} else if (species == "mouse") {
  newNames$speciesOneToOne <- ifelse(is.na(newNames$mouseOneToOne), newNames$uniqueHier, newNames$mouseOneToOne)
} else if (species == "woodchuck") {
  newNames$speciesOneToOne <- newNames$uniqueHier
}
# Grab dataframe
woodchuckClusterAverages <- woodchuckClusterAverages$SCT
# Replace names with one-to-one orthologue of particular species
row.names(woodchuckClusterAverages) <- newNames$speciesOneToOne
# Make sure formatted correctly
woodchuckClusterAverages <- as.data.frame(woodchuckClusterAverages)
# Order by gene name
woodchuckClusterAverages <- woodchuckClusterAverages[order(row.names(woodchuckClusterAverages)),]
# Sanity check
head(woodchuckClusterAverages)

Dataset-specific setup

Correlation of woodchuck PBMCs with human 68k PBMC dataset from 10X Genomics

# Sanity check
head(allCellsMatrix)

       Activated CD8+ Naive CD8+ Memory and Reg T Naive CD4+          NK       CD8+          B
42430      -0.7290166 -1.4467216      -0.19073790 -0.4832807  0.02379348 -0.8928406 -0.2258430
42431      -0.2300094 -0.7300300      -0.26479348 -0.7995981 -0.18652940 -0.3474056 -0.6952459
42618      -0.1212209 -0.6309294      -0.09745206 -0.4354971 -0.36419075 -0.8527715 -0.3826776
A4GALT     -0.4058024 -1.0280327      -0.68648235 -0.4734361 -0.11159565 -0.5478332 -0.1352675
AATK       -1.0033368 -0.9033733      -0.11847519 -1.0810861 -0.10366579 -0.4257702 -0.2369504
ABCA1      -0.6939429 -1.0596759      -0.79368826 -0.8411861 -0.88868386 -0.1144697  0.6312458
       Megakaryocytes Monocytes and Dendritic B, Dendritic, T
42430        1.314882              1.50210970       1.1276549
42431        2.578802              0.03087083       0.6439395
42618        2.691420             -0.14763062       0.3409502
A4GALT       2.062829             -0.23671805       1.5623391
AATK         1.958543              1.15513309       0.7589817
ABCA1        1.186970              1.02547753       1.5479533

Correlation of woodchuck liver with human liver dataset from MacParland et al. (2018)

Correlation of woodchuck liver with human liver dataset from Aizarani et al.

allCellsMatrix <- allCellsMatrix[,as.character(sort(as.numeric(colnames(allCellsMatrix))))]
# Rename columns to be more meaningful (not totally confident I got all correct)
colnames(allCellsMatrix) <- c('NK, NKT and T cells (1)',
                              'Kupffer cells (2)',
                              'NK, NKT and T cells (3)',
                              'EPCAM+ cells and cholangiocytes (4)',
                              'NK, NKT and T cells (5)',
                              'Kupffer cells (6)',
                              'EPCAM+ cells and cholangiocytes (7)',
                              'B cells (8)',
                              'Liver sinusoidal endothelial cells (9)',
                              'Macrovascular endothelial cells (10)',
                              'Hepatocytes (11)',
                              'NK, NKT and T cells (12)',
                              'Liver sinusoidal endothelial cells (13)',
                              'Hepatocytes (14)',
                              'Other endothelial cells (15)',
                              'Other (16)',
                              'Hepatocytes (17)',
                              'NK, NKT and T cells (18)',
                              'NK, NKT and T cells (19)',
                              'Liver sinusoidal endothelial cells (20)',
                              'Macrovascular endothelial cells (21)',
                              'B cells (22)',
                              'Kupffer cells (23)',
                              'EPCAM+ cells and cholangiocytes (24)',
                              'Kupffer cells (25)',
                              'Other endothelial cells (26)',
                              'Other (27)',
                              'NK, NKT and T cells (28)',
                              'Macrovascular endothelial cells (29)',
                              'Hepatocytes (30)',
                              'Kupffer cells (31)',
                              'Macrovascular endothelial cells (32)',
                              'Stellate cells and myofibroblasts (33)',
                              'B cells (34)',
                              'Other endothelial cells (35)',
                              'Other (36)',
                              'Other (37)',
                              'B cells (38)',
                              'EPCAM+ cells and cholangiocytes (39)')
# Order by row name
allCellsMatrix <- allCellsMatrix[order(row.names(allCellsMatrix)),]
# Sanity check
head(allCellsMatrix)

speciesData <- "aizarani"

Correlation of woodchuck liver with woodchuck PBMCs. For this correlation, read in the woodchuck liver dataset at the beginning of this script and then read in the woodchuck PBMCs below

groups <- "general_cell_labels"
# Start with liver and read in woodchuck PBMCs again
load("~/Dropbox/Zoe/scf_version/analysis/healthy_sc/seurat_objects/no_dropletQC/integrated_PBMC_cca_kanchor5_scClustViz.RData")
Idents(scSeurat) <- "integrated_snn_res.0.6"
# Find cluster averages
pbmcClusterAverages <- AverageExpression(scSeurat,
                                         assays = "SCT",
                                         slot = "scale.data",
                                         group.by = groups)
pbmcClusterAverages <- as.data.frame(na.omit(t(scale(t(as.matrix(pbmcClusterAverages$SCT))))))
# Order by row name
allCellsMatrix <- pbmcClusterAverages[order(row.names(pbmcClusterAverages)),]
speciesData <- "PBMC"

Output plots for all correlations

# Now find intersecting genes
matches <- intersect(row.names(allCellsMatrix),
                     row.names(woodchuckClusterAverages))
# Look at how many genes matched
length(matches)

[1] 165

# Make new matrices with only matching gene names
toCor <- allCellsMatrix[matches,]
woodchuckAveragesCor <- woodchuckClusterAverages[matches,]
# Do Pearson
pearVal <- cor(toCor, woodchuckAveragesCor, method = "pearson")
heatmap(pearVal)
        #main = paste("Pearson correlation of", speciesData, "vs", woodchuckData),
        #xlab = woodchuckData,
        #ylab = speciesData)
        #margins = c(6,11))
#Rowv = NA,
#Colv = NA)
pdf(paste("./figures/", tissue, "/", tissue, "_", speciesData, "_PearsonCor.pdf", sep = ""),
    height = 13, width = 13)

heatmap(pearVal, margins = c(13,13))
dev.off()

png 
  2 

# Do Spearman
spearVal <- cor(toCor, woodchuckAveragesCor, method = "spearman")
heatmap(spearVal)
        #main = paste("Spearman correlation of", speciesData, "vs", woodchuckData),
        #xlab = woodchuckData,
        #ylab = speciesData)
        #margins = c(6,11))
#Rowv = NA,
#Colv = NA)
pdf(paste("./figures/", tissue, "/", tissue, "_", speciesData, "_SpearmanCor.pdf", sep = ""),
    height = 13, width = 13)

heatmap(spearVal, margins = c(13,13))
dev.off()

png 
  2